coef.names = c('Intercept', 'Segregation', 'Outgroup Proportion',  'Segregation x Outgroup Proportion', 'Male', 'Age', 'Mixed Ethnicity', 'Other Ethnicity', 'Sephardic Ethnicity', 'Political Ideology','High Income','Low Income','College','Immigrant','Arab Population','Jerusalem')

##cycle through different dats, with different selection controls, all of these are people that are likley not to be selection problems: 1) people that do not prefer to live exclusively with the ingroup, if they had the choice, 2) people that are not satisified with where they live, 3) people that do not feel they are free to move, 4) low income people, 5) people who do not list religious people as a reason affecting city choice, 6) people who do not list religious people as a reason affecting neighborhood choice
dats.uo = list(
  dat.uo[dat.uo$ingroup.neighborhood.prefer.pref == 0,], 
	dat.uo[dat.uo$satisfaction.index.high == 0,],
	dat.uo[dat.uo$free.to.move == 0,],
	dat.uo[dat.uo$low.income == 1,],
	dat.uo[dat.uo$move.religious.city == 0,],
	dat.uo[dat.uo$move.religious.nbhood == 0,],
	dat.uo[dat.uo$selecters == 0,],
  dat.uo[dat.uo$outgroup.houses.exclude == 0,],
  dat.uo ##end with complete sample
  )
  
  dats.non.uo = list(
  dat.non.uo[dat.non.uo$ingroup.neighborhood.prefer.pref == 0,],
	dat.non.uo[dat.non.uo$satisfaction.index.high == 0,],
	dat.non.uo[dat.non.uo$free.to.move == 0,],
	dat.non.uo[dat.non.uo$low.income == 1,],
	dat.non.uo[dat.non.uo$move.religious.city == 0,],
	dat.non.uo[dat.non.uo$move.religious.nbhood == 0,],
	dat.non.uo[dat.non.uo$selecters == 0,],
  dat.non.uo[dat.non.uo$outgroup.houses.exclude == 0,],  
  dat.non.uo ##end with complete sample
  )

##lists to store output
uo.dict.list = list()
sec.dict.list = list()

uo.pub.list = list()
sec.pub.list = list()

uo.task.list = list()
sec.task.list = list()


for(i in 1:length(dats.uo)){
	master.dat.uo = dats.uo[[i]]
	master.dat.non.uo = dats.non.uo[[i]]
			   
	##############################################		
	###need complete cases so subset below		

  master.use.vars = c('outgroup.dissim.yeshiva','outgroup.yeshiva','demo1.sex','demo1.age','demo1.ethnicity','demo2.left_right','index.city','nonjewish_pcnt','income','college','immigrant','jerusalem')

		   
	##master formula, combine with the dependent variable in each script
  master.formula = 'outgroup.dissim.yeshiva*outgroup.yeshiva+demo1.sex+demo1.age+demo1.ethnicity+demo2.left_right+income+college+immigrant+nonjewish_pcnt+jerusalem'


	if(i == 4){ ##drop income from formula when subsetting to low income
	  master.formula ='outgroup.dissim.yeshiva*outgroup.yeshiva+demo1.sex+demo1.age+demo1.ethnicity+demo2.left_right+college+immigrant+nonjewish_pcnt+jerusalem'
	  	}

  ##UO
	use.vars = c('in.out.diff.dictator',master.use.vars)
	use.dat = master.dat.uo[,use.vars]
	use.dat = use.dat[complete.cases(use.dat),]

	reg.uo.dict = lm(as.formula(paste('in.out.diff.dictator~',master.formula,sep='')), data = use.dat)
	clustered.se = cl(use.dat, reg.uo.dict, use.dat$index.city)
	reg.uo.dict$se = clustered.se
	uo.dict.list[[i]] = reg.uo.dict

	
	use.vars = c('public.ingroup.preference',master.use.vars)
	use.dat = master.dat.uo[,use.vars]
	use.dat = use.dat[complete.cases(use.dat),]

	reg.uo.pub = lm(as.formula(paste('public.ingroup.preference~',master.formula,sep='')), data = use.dat)
	clustered.se = cl(use.dat, reg.uo.pub, use.dat$index.city)
	reg.uo.pub$se = clustered.se
	uo.pub.list[[i]] = reg.uo.pub
	
  
	use.vars = c('ingroup.task',master.use.vars)
	use.dat = master.dat.uo[,use.vars]
	use.dat = use.dat[complete.cases(use.dat),]

	reg.uo.task = lm(as.formula(paste('ingroup.task~',master.formula,sep='')), data = use.dat)
	clustered.se = cl(use.dat, reg.uo.task, use.dat$index.city)
  reg.uo.task$se = clustered.se
	uo.task.list[[i]] = reg.uo.task
	
  
  ##Non UO
	use.vars = c('in.out.diff.dictator',master.use.vars)
	use.dat = master.dat.non.uo[,use.vars]
	use.dat = use.dat[complete.cases(use.dat),]

	reg.sec.dict = lm(as.formula(paste('in.out.diff.dictator~',master.formula,sep='')), data = use.dat)
	clustered.se = cl(use.dat, reg.sec.dict, use.dat$index.city)
	reg.sec.dict$se = clustered.se
	sec.dict.list[[i]] = reg.sec.dict

	
	use.vars = c('public.ingroup.preference',master.use.vars)
	use.dat = master.dat.non.uo[,use.vars]
	use.dat = use.dat[complete.cases(use.dat),]

	reg.sec.pub = lm(as.formula(paste('public.ingroup.preference~',master.formula,sep='')), data = use.dat)
	clustered.se = cl(use.dat, reg.sec.pub, use.dat$index.city)
	reg.sec.pub$se = clustered.se
  sec.pub.list[[i]] = reg.sec.pub
  
	use.vars = c('ingroup.task',master.use.vars)
	use.dat = master.dat.non.uo[,use.vars]
	use.dat = use.dat[complete.cases(use.dat),]

	reg.sec.task = lm(as.formula(paste('ingroup.task~',master.formula,sep='')), data = use.dat)
	clustered.se = cl(use.dat, reg.sec.task, use.dat$index.city)
	reg.sec.task$se = clustered.se
  sec.task.list[[i]] = reg.sec.task
  
	}


##write out tables:
##dictator
outtable = apsrtable(uo.dict.list[[1]], uo.dict.list[[2]], uo.dict.list[[3]], uo.dict.list[[4]], uo.dict.list[[5]], uo.dict.list[[6]], uo.dict.list[[7]], uo.dict.list[[8]],
		Sweave = T,
		coef.names = coef.names,
	notes = ''
	)
writeLines(
  outtable, 'output/appendix/Table_A9.tex', sep = '')

outtable = apsrtable(sec.dict.list[[1]], sec.dict.list[[2]], sec.dict.list[[3]], sec.dict.list[[4]], sec.dict.list[[5]], sec.dict.list[[6]], sec.dict.list[[7]], sec.dict.list[[8]],
		Sweave = T,
		coef.names = coef.names,
	notes = ''
	)
writeLines(
	outtable, 'output/appendix/Table_A10.tex', sep = '')


##task
outtable = apsrtable(uo.task.list[[1]], uo.task.list[[2]], uo.task.list[[3]], uo.task.list[[4]], uo.task.list[[5]], uo.task.list[[6]], uo.task.list[[7]], uo.task.list[[8]],
                     Sweave = T,
                     coef.names = coef.names,
                     notes = ''
                        )
writeLines(
  outtable, 'output/appendix/Table_A11.tex', sep = '')

outtable = apsrtable(sec.task.list[[1]], sec.task.list[[2]], sec.task.list[[3]], sec.task.list[[4]], sec.task.list[[5]], sec.task.list[[6]], sec.task.list[[7]], sec.task.list[[8]],
                     Sweave = T,
                     coef.names = coef.names,
                     notes = ''
                      )
writeLines(
  outtable, 'output/appendix/Table_A12.tex', sep = '')


##public
outtable = apsrtable(uo.pub.list[[1]], uo.pub.list[[2]], uo.pub.list[[3]], uo.pub.list[[4]], uo.pub.list[[5]], uo.pub.list[[6]], uo.pub.list[[7]], uo.pub.list[[8]],
                     Sweave = T,
                     coef.names = coef.names,
                     notes = ''
                    )
writeLines(
    outtable, 'output/appendix/Table_A13.tex', sep = '')
  
outtable = apsrtable(sec.pub.list[[1]], sec.pub.list[[2]], sec.pub.list[[3]], sec.pub.list[[4]], sec.pub.list[[5]], sec.pub.list[[6]], sec.pub.list[[7]], sec.pub.list[[8]],
                     Sweave = T,
                     coef.names = coef.names,
                     notes = ''
                  )
writeLines(
    outtable, 'output/appendix/Table_A14.tex', sep = '')




######################graph coefficients
##tables for storing coefficients and se's
coef.table.outgroup = matrix(ncol = 4, nrow = length(dats.uo))
coef.table.dissim = matrix(ncol = 4, nrow = length(dats.uo))
coef.table.interaction = matrix(ncol = 4, nrow = length(dats.uo))
N.table = matrix(ncol = 2, nrow = length(dats.uo)) ##store length

for(i in 1:length(dats.uo)){  ##extract coefficients
  coef.table.outgroup[i,1] = summary(sec.dict.list[[i]])$coefficients['outgroup.dissim.yeshiva','Estimate']
  coef.table.outgroup[i,2] = summary(sec.dict.list[[i]])$coefficients['outgroup.dissim.yeshiva','Std. Error']
  coef.table.outgroup[i,3] = summary(uo.dict.list[[i]])$coefficients['outgroup.dissim.yeshiva','Estimate']
  coef.table.outgroup[i,4] = summary(uo.dict.list[[i]])$coefficients['outgroup.dissim.yeshiva','Std. Error']
    
  coef.table.dissim[i,1] = summary(sec.dict.list[[i]])$coefficients['outgroup.yeshiva','Estimate']
  coef.table.dissim[i,2] = summary(sec.dict.list[[i]])$coefficients['outgroup.yeshiva','Std. Error']
  coef.table.dissim[i,3] = summary(uo.dict.list[[i]])$coefficients['outgroup.yeshiva','Estimate']
  coef.table.dissim[i,4] = summary(uo.dict.list[[i]])$coefficients['outgroup.yeshiva','Std. Error']
  
  coef.table.interaction[i,1] = summary(sec.dict.list[[i]])$coefficients['outgroup.dissim.yeshiva:outgroup.yeshiva','Estimate']
  coef.table.interaction[i,2] = summary(sec.dict.list[[i]])$coefficients['outgroup.dissim.yeshiva:outgroup.yeshiva','Std. Error']
  coef.table.interaction[i,3] = summary(uo.dict.list[[i]])$coefficients['outgroup.dissim.yeshiva:outgroup.yeshiva','Estimate']
  coef.table.interaction[i,4] = summary(uo.dict.list[[i]])$coefficients['outgroup.dissim.yeshiva:outgroup.yeshiva','Std. Error']
  
  N.table[i,1]  = nrow(dats.uo[[i]])
  N.table[i,2]  = nrow(dats.non.uo[[i]])
  
  }


##make plots, one for UO, one for secular

##graphical parameters
ylims = c(.75,length(dats.uo)+.5)
outgroup.ys = 1:length(dats.uo)
dissim.ys = outgroup.ys + .2
interaction.ys = outgroup.ys + .4
text.ys = outgroup.ys + .6
point.cex = 1
line.lwd = 1
line.offset = .025
lab.cex = .5
text.cex = .55


##labels for text
text.labels = c('do not prefer homogenous neighborhood', 'not satisfied with location', 'not free to move', 'low income',
                  'religion not a factor in city choice','religion not a factor in neighborhood choice', 'did not select into high-ingroup area',
                'ideal neighborhood not homogenous', 'all respondents')


##sec specific
##confident interval mats
cis.outgroup = coef.table.outgroup[,2]*1.96
cis.dissim = coef.table.dissim[,2]*1.96
cis.interaction = coef.table.interaction[,2]*1.96
sec.xlims = c(-60,200)
sec.labels = paste(text.labels, " (",N.table[,1],")",sep = "") ##paste N into labels for chart

pdf('output/Figure_4a.pdf',
      width = 3.25,
      height = 7.5)  
par(bty = 'n')
par(mar = c(4,1,1,1))
plot(0,0,
     xlim = sec.xlims,
     ylim = ylims,
     type = 'n',
     yaxt = 'n',
     ylab = '',
     xlab = 'Regression Coefficients',
     cex.lab = lab.cex,
     cex.axis = lab.cex
  )


##make vertical lines to represent the points of the main regressions
abline(v = coef.table.outgroup[length(dats.uo),1],
        col = 'gray',
       lty = 2,
       lwd = line.lwd+.5)
abline(v = coef.table.dissim[length(dats.uo),1],
       col = 'gray',
       lty = 2,
       lwd = line.lwd+.5)
abline(v = coef.table.interaction[length(dats.uo),1],
       col = 'gray',
       lty = 2,
       lwd = line.lwd+.5)


points(coef.table.outgroup[,1],
      outgroup.ys,
       pch = 15,
       cex = point.cex,
       col = line.color)
points(coef.table.dissim[,1],
       dissim.ys,
       pch = 16,
       cex = point.cex,
       col = line.color)
points(coef.table.interaction[,1],
       interaction.ys,
       pch = 17,
       cex = point.cex,
       col = line.color)

##confidence intervals
for(i in 1:length(cis.outgroup)){
   lines(x = c(coef.table.outgroup[i,1] - cis.outgroup[i],coef.table.outgroup[i,1] + cis.outgroup[i]), 
         y = c(outgroup.ys[i],outgroup.ys[i]),
         lty = 2,
         lwd = line.lwd,
         col = line.color)
   lines(x = c(coef.table.dissim[i,1] - cis.dissim[i],coef.table.dissim[i,1] + cis.dissim[i]), 
         y = c(dissim.ys[i],dissim.ys[i]),
         lty = 2,
         lwd = line.lwd,
         col = line.color)
   lines(x = c(coef.table.interaction[i,1] - cis.interaction[i],coef.table.interaction[i,1] + cis.interaction[i]), 
         y = c(interaction.ys[i],interaction.ys[i]),
         lty = 2,
         lwd = line.lwd,
         col = line.color)  
     }
 text(sum(sec.xlims)/2,
      text.ys,
      sec.labels,
      cex = text.cex,
      col = line.color)
dev.off()



##UO specific
##confident interval mats
cis.outgroup = coef.table.outgroup[,4]*1.96
cis.dissim = coef.table.dissim[,4]*1.96
cis.interaction = coef.table.interaction[,4]*1.96
uo.xlims = c(-140,140)
uo.labels = paste(text.labels, " (",N.table[,2],")",sep = "") ##paste N into labels for chart

pdf('output/Figure_4b.pdf',
    width = 3.25,
    height = 7.5)  
par(bty = 'n')
par(mar = c(4,1,1,1))
plot(0,0,
     xlim = uo.xlims,
     ylim = ylims,
     type = 'n',
     yaxt = 'n',
     ylab = '',
     xlab = 'Regression Coefficients',
     cex.lab = lab.cex,
     cex.axis = lab.cex
  )

##make vertical lines to represent the points of the main regressions
abline(v = coef.table.outgroup[length(dats.uo),3],
       col = 'gray',
       lty = 2,
       lwd = line.lwd+.5)
abline(v = coef.table.dissim[length(dats.uo),3],
       col = 'gray',
       lty = 2,
       lwd = line.lwd+.5)
abline(v = coef.table.interaction[length(dats.uo),3],
       col = 'gray',
       lty = 2,
       lwd = line.lwd+.5)


points(coef.table.outgroup[,3],
       outgroup.ys,
       pch = 15,
       cex = point.cex,
       col = line.color)
points(coef.table.dissim[,3],
       dissim.ys,
       pch = 16,
       cex = point.cex,
       col = line.color)
points(coef.table.interaction[,3],
       interaction.ys,
       pch = 17,
       cex = point.cex,
       col = line.color)

##confidence intervals
for(i in 1:length(cis.outgroup)){
  lines(x = c(coef.table.outgroup[i,3] - cis.outgroup[i],coef.table.outgroup[i,3] + cis.outgroup[i]), 
        y = c(outgroup.ys[i],outgroup.ys[i]),
        lty = 2,
        lwd = line.lwd,
        col = line.color)
  lines(x = c(coef.table.dissim[i,3] - cis.dissim[i],coef.table.dissim[i,3] + cis.dissim[i]), 
        y = c(dissim.ys[i],dissim.ys[i]),
        lty = 2,
        lwd = line.lwd,
        col = line.color)
  lines(x = c(coef.table.interaction[i,3] - cis.interaction[i],coef.table.interaction[i,3] + cis.interaction[i]), 
        y = c(interaction.ys[i],interaction.ys[i]),
        lty = 2,
        lwd = line.lwd,
        col = line.color)  
}
text(sum(uo.xlims)/2,
     text.ys,
     uo.labels,
     cex = text.cex,
     col = line.color)
dev.off()
